Method for time lapse reservoir monitoring

ABSTRACT

A method of comparing multiple seismic survey data sets of a reservoir, is provided, wherein a first seismic survey data set is taken at a first time and a second seismic survey data set is taken at a second time, to detect changes in the reservoir between the first time and the second time. Generally, the method comprises design of a crossequalization function for the second data set based on a comparison of unchanging portions of the two data sets. Also provided is a processing method for preparing each survey according to similar processing steps with information taken from each survey.

This application is a continuation-in-part of U.S. application Ser. No.08/713,948, Sep. 13, 1996, abn, assigned to the assignee of the presentinvention, and including an inventor common to the inventors of thepresent invention.

BACKGROUND OF THE INVENTION

This invention relates to oil and gas reservoir management and, morespecifically, to time lapse reservoir seismic signal processing.

Reservoir characterization and monitoring in the oil and gas field areimportant parts of reservoir management and hydrocarbon production.Effective reservoir management is a major goal of energy producingcompanies as they try to reduce finding costs, optimize drillinglocations and increase financial returns. One technique that isattempted in this endeavor is time-lapse (also known as 4D) seismicmonitoring. As fluids are extracted, swept, or injected throughproduction and recovery, changes in the effective elastic properties ofthe reservoir rocks occur. The ability to monitor reservoir changes as afunction of time by the use of seismic methods can lead to betterlocation of production and infill wells, the possibility of locatingunswept zones, and more efficient field maintenance, thus raising theoverall value of the production lease.

In a two-dimensional approach, seismic monitoring has been examined incrosswell procedures. However, the repeated results have only beenqualitatively compared, and two dimensional time-lapse surveys, so far,do not contain the type of information desired in modern reservoirmanagement (for example, see Paulsson, et. al, 1994, The Leading Edge,incorporated herein by reference). Time-lapse 3D techniques have alsobeen attempted, but they involve complex modeling procedures and requirea great deal of processing without using direct data available in thesurvey itself. See, e.g., U.S. Pat. No. 4,969,130, incorporated hereinby reference.

One problem in time lapse processing is that many conditions change overtime, not just the changes in the state of the reservoir. For example,the locations of the source and receiver in the second survey willnecessarily be different from those in the first. Further, the tide in asecond marine survey may be higher or lower, as may the temperature ofthe air and water. Likewise, the specific characteristics of the sourcesand receivers used in the second survey will be different. Otherdifferences, besides changes in reservoir state also occur, such asdifferences in the manner in which the two surveys are processed. Thus,there is a need for a method of dealing with the two surveys wherebyprocessing differences do not detrimentally affect the result of thecomparison.

For example, in gathering seismic data, a source is used to generateseismic waves which reflect from the reflectors in the earth (e.g. layerboundaries) and are received at receivers. In some cases, the sourcesignature is a spike, although, in reality, it is not perfect. Duringits journey through the strata and reflectors, the signal shape ischanged, and the reflection signal received at the receivers is,therefore, no longer a spike, or even close. Deconvolution is theprocess by which the shape of the reflection signal is “whitened” torecreate the spike shape of the data.

In another example, a broad, band-limited signal is used, which is zerophase. Deconvolution is used in such a case to remove the distortionscaused by the earth.

In still another example, in performing the deconvolution in thefrequency domain, all frequency samples are multiplied to bring them toan equal level, following the assumption that the source is a minimumphase signal, immediately rising to a peak and then dying. This isaccomplished by autocorrelating the trace in the time domain multipletimes, at a series of lag samples, which results in a generallysymmetrical wavelet. The power spectrum of the wavelet is then analyzed,to determine what multipliers are needed at each frequency to flattenthe frequency spectrum. This process is performed on a windowed basis,both along each trace, and across the record (as used herein, the term“record” refers to, alternatively for example, a common receiver record,a CMP record, a common shot record, a stacked trace record, etc.) Theautocorrelation is performed on various windows, and the results areaveraged to give the spectrum. From that spectrum, the operators neededto flatten the spectrum are chosen. The operators are then applied toall of the traces used in the input. Typically, the window is about 10times the length of the operator to be generated, measured in number ofsamples. The deconvolution process and the design of a deconvolutionoperator are well known in the art, and it is not limited to thefrequency domain example, above. It is also routinely performed in thetime domain. See, e.g., Yilmaz, Investigations in Geophysics, Vol. 2,Seismic Data Processing, Society of Exploration Geophysicists (1987) andreferences cited therein.

In performing deconvolution, it is important to design a deconvolutionoperator dependent upon the data of the survey, in order to account forthe specific source signature, and other equipment distortions thatoccur in the data. Therefore, the data in each survey has been adjustedthrough the use of a specific optimum deconvolution operator which isnot applicable to other surveys. The result of this difference in theuse of the separate deconvolution operators in time-lapse surveys isthat structure appears in the difference records when two records aresubtracted. This result is undesirable. However, to date, no one hasproposed a practical solution to the problem.

SUMMARY OF THE INVENTION

It is an object of the present invention to address the above problems.

It has been found that, contrary to earlier beliefs, a singledeconvolution operator can be used on multiple sets of data, not onlywithout detrimental results, but improving the quality of the processingof time-lapse comparisons of seismic surveys. Accordingly, in one aspectof the present invention, a method of deconvolution of multiple sets ofseismic data from the same geographic area is provided, the methodcomprising: designing of a deconvolution operator dependent upon datafrom at least two of the sets of seismic data, wherein the at least twosets of seismic data were recorded at different times or calendar dates;applying the deconvolution operator in a deconvolution process to bothof the at least two sets of data; and conducting further time-lapseprocessing to form a difference record.

According to one embodiment of the invention, the conducting of furthertime-lapse processing comprises: providing a first reflection event (forexample, a wavelet) in the first seismic survey data set having acorresponding second reflection event in the second seismic survey dataset, wherein the first reflection event and the second reflection eventrepresent an unchanged portion of geologic structure in or near thereservoir and wherein the first reflection event is represented by afirst set of event parameters and the second reflection event isrepresented by a second set of event parameters. Next, an acceptancethreshold difference function between the first set of event parametersand the second set of event parameters is provided. Then, acrossequalization function is determined to apply to the second set ofevent parameters.

According to another aspect of the invention, the crossequalizationfunction is determined such that, upon application of thecrossequalization function to the second set of event parameters, acrossequalized set of event parameters is defined, and the differencebetween the first set of event parameters and the crossequalized set ofevent parameters is below the threshold difference function. Next, thecrossequalization function is applied to a third reflection event, thethird reflection event being related to the second data set, wherein acrossequalized third reflection event is defined, wherein the thirdreflection event has a corresponding fourth reflection event in thefirst data set, and wherein the third and fourth reflection eventsrepresent a changing portion of the reservoir.

Comparison of the crossequalized third reflection event to the fourthreflection event by subtracting the crossequalized third reflectionevent from the fourth reflection event results in the desiredinformation.

According to a more specific example embodiment, said providing saidacceptance threshold difference function comprises: iterative selectionof event parameter modifications to the second set of event parameters,application of the event parameter modifications to the second set ofevent parameters, wherein a modified set of event parameters is defined,comparison of the modified set of event parameters to the first set ofevent parameters, wherein said iterative selection continues until aconvergence is reached, and wherein the acceptance threshold differencefunction comprises the modified set of event parameters at convergence.Example event parameters comprise any combination of amplitude, phase,bandwidth, and time, or any of the foregoing individually.

According to another example embodiment, the determining of thecrossequalization function comprises: iterative selection of eventparameter modifications to the second set of event parameters,application of the event parameter modifications to the second set ofevent parameters, wherein a modified set of event parameters is defined,comparison of the modified set of event parameters to the first set ofevent parameters, and providing an acceptance threshold difference,wherein said iterative selection continues until a comparison resultfrom said comparison designates a difference between the first set ofevent parameters and the modified set of event parameters below theacceptance threshold difference.

According to yet another example embodiment, said providing anacceptance threshold difference function comprises: providing a windowedtrace difference between a time window of a first trace from the firstseismic survey data set and a time window of a second trace from thesecond seismic survey data set, wherein the second trace includesreflection events corresponding to reflection events in the first traceand wherein the time window of the second trace is substantially thesame as the time window of the first trace, and providing a ratio of thewindowed trace difference over the time window of the first trace, andchoosing the acceptance threshold difference to be less than the ratio.

The time windows in both the unchanging and changing portions of thereservoir have similar spectral characteristics. For example, if thedata from the reservoir has a dominant frequency of 30 Hz, the timewindow used should be picked from an unchanging portion of the surveyhaving a dominant frequency as close to 30 Hz as possible. Likewise,phase changes in the reservoir and the unchanging portion should be asclose as possible. It is preferred, however, to err on the side ofbroader time windows. For example, if the reservoir dominant frequencyis 30 Hz, a window having 35 Hz is considered preferable to one of 25Hz. Such bandwidth error of less than about 25% in frequency bandwidthis believed to yield adequate results. Best results should be seen whenbandwidth error is below 10%.

In still another embodiment, said providing an acceptance thresholddifference function comprises: providing a windowed trace differencebetween a time window of the square of a first trace from the firstseismic survey data set and a time window of the square of a secondtrace from the second seismic survey data set, wherein the second traceincludes reflection events corresponding to reflection events in thefirst trace and wherein the time window of the second trace issubstantially the same as the time window of the first trace, andproviding a ratio of the windowed trace difference over the time windowof the square of the first trace, choosing the acceptance thresholddifference to be less than the ratio.

In still another embodiment, said applying the crossequalizationfunction to a third reflection event in the second data set comprisesconvolution between the crossequalization function and the thirdreflection event in the second data set, said first data set comprises atrace from a seismic receiver. Alternatively, said first data set andsaid second data set comprise a summed set of traces from a set ofseismic receivers, or CMP (“common mid-point”) data. In still a furtherembodiment, said first data set and said second data set comprise shotdata. In further alternatives, said first data set and said second dataset comprise prestack data or migrated data.

In many embodiments, said first data set and said second data set aresubjected to equivalent prestack processes. For example, in addition tothe deconvolution described above, in some embodiments the first dataset uses the same designature process as the second data set, the samenoise attenuation processing steps as second data set, and the samemultiple attenuation processing as the second data set. Further, in manyembodiments, the same DMO operator is used on first data set as on thesecond data set, and migration on the first data set is conducted withthe same velocity field as migration on the second data set.

Finally according to a further aspect of the invention, a method isprovided for performing time-lapse seismic survey signal processing, themethod comprising: performing a set of processing steps on the firstsurvey; performing the set of processing steps on the second survey,wherein the set of processing steps is dependent upon a set of seismicsignal parameters; choosing at least one of the set of parameters by aselection process dependent upon data from both surveys; and (b)applying the at least one of the sets of parameters in the at least oneof the set of processing steps to both the first survey and the secondsurvey.

DESCRIPTION OF THE DRAWINGS

For a more complete understanding of the present invention and forfurther advantages thereof, reference is made to the followingDescription of Example Embodiments of the Invention, taken inconjunction with the accompanying drawings, in which:

FIG. 1A is a plot of a first record of a survey taken at a first time,without deconvolution.

FIG. 1B is a plot of the first record of FIG. 1A with deconvolution,wherein the deconvolution operator was designed dependent upon the dataof the first record.

FIG. 1C is a plot of the raw data of FIG. 1A with deconvolution applied,wherein the deconvolution operator was designed dependent on both thedata of FIG. 1A and another survey.

FIG. 2A is a plot of a second record of a second survey taken of thesame geographic area as that of FIG. 1A, but at a different time,without deconvolution.

FIG. 2B is a plot of the second record of FIG. 2A, with deconvolution,wherein the deconvolution operator was designed dependent upon the dataof the second record.

FIG. 2C is a plot of the second record of FIG. 2A, with deconvolution,wherein the deconvolution operator was designed dependent upon the dataof both the record of FIG. 1A and the record of FIG. 2A.

FIG. 3A is a plot of the difference between FIGS. 1A and 2C.

FIG. 3B is a plot of the difference between FIGS. 1B and 2B

FIG. 3C is a plot of the difference between FIGS. 1C and 2C

It is to be noted, however, that the appended drawings illustrate onlytypical embodiments of this invention and are therefore not to beconsidered limiting of its scope, for the invention may admit to otherequally effective embodiments.

DESCRIPTION OF EXAMPLE EMBODIMENTS OF THE INVENTION

As used herein, “crossequalization” is the catch-all industry term forthe match-filtering, amplitude scaling and static corrections necessaryfor time-lapse seismic monitoring. In essence, a wavelet operator oroperators are estimated to shape and match the reflection data from onesurvey to another. Typically, operators are designed over staticreflectors that exclude the reservoir where meaningful changes in (pore)fluid state may occur. Theoretically, the difference between two volumesof data after crossequalization should be zero everywhere (inside andoutside the design window) except where there are changes in thereservoir. All of the static or non-reservoir events should besubtracted out, leaving only changes in dynamic events (i.e., fluidreplacement within the pores).

Whether a change in the reservoir state is seismically detectable isdetermined by the reservoir's lithology and the acoustic properties ofthe pore fluids (which are temperature and pressure dependent), the typeof production or recovery process, the resolution (spatial and temporal)of the data, in addition to the repeatability of the seismic system. Incertain reservoirs with gas drives, gas injection may reduce theacoustic impedance (oil being displaced by gas) sufficiently to induce abright spot. See, e.g., Fulp and Greaves, GEOPHYSICS, 1987, incorporatedherein by reference, for a bright spot generated by a fire flood. Inothers, gas injection or water injection (water drive) may have theopposite reaction and generate a dim spot. Where steam flooding isinvolved, velocity pull-down or sag often is used to indicate steamfronts to determine to what extent the recovery has progressed. Forexample, see Lumley, et. al., SEG Expanded Abstracts, 1995, incorporatedherein by reference. As aforementioned, petrophysical and seismicinterpretation are required to understand what attribute is meaningful,and how much change one might observe.

Processing

The rudimentary data processing steps toward the goal of obtainingspatially-correct, relative amplitude preserved data include NMO/DMO andprestack (zero-offset) migration used for the forthcoming examples.Other appropriate prestack procedures will occur to those of skill inthe art, including data compression and sampling techniques that can beemployed to reduce the computational overhead and preserve the qualityof the prestack processing that will occur to those of skill in the art.After such processing, relative amplitude data conditioning yieldsprestack migrated CMP data on which crossequalization and time-lapsedifferencing measures are performed. Depending on the producingreservoir, according to alternative embodiments of the invention, CMPs,AVO attributes, and/or migrated stacks from each survey arecrossequalized and subtracted (differenced) to observe interpore fluiddisplacements.

It is recommended that each survey be processed as close to identicallyas possible (which is not normally done in time-lapse processing).Before the present invention, whatever processed data existed onworkstations (with dissimilar preconditioning procedures and algorithms)was used, and removal of the resultant variations in data was attemptedthrough crossequalization. Better crossequalizations and more meaningfulseismic differences will be obtained by reprocessing existing data withthe same methodology and software as the recent seismic monitoringsurvey. Thus, it is particularly preferred to designature, deconvolve,attenuate noise and multiples with the same parameterization, survey tosurvey, just as it is important to use the same DMO operator and migratethe data with the same velocity field, to obtain premier time-lapsemeasurements.

Therefore, according to one aspect of the invention, a method isprovided for performing time-lapse seismic survey signal processing,wherein a set of processing steps is performed on the first survey, theset of processing steps in performed on the second survey, and the setof processing steps is dependent upon a set of seismic signalparameters, the method comprising: (a) choosing at least one of the setof parameters by a selection process dependent upon data from bothsurveys; and (b) applying the at least one of the sets of parameters inthe at least one of the set of processing steps to both the first surveyand the second survey.

For example, it has been found that improved deconvolution processing ofthe two surveys results from the design of a deconvolution operatordependent upon multiple sets of data. This deconvolution operator isthen applied to both data sets. The result of such a design, whileperhaps not optimum for either particular set, is, nevertheless, optimumfor time-lapse comparison purposes.

According to one example embodiment of this aspect of the invention, amethod of deconvolution of multiple sets of seismic data from the samegeographic area is provided, the method comprising: designing of adeconvolution operator dependent upon data from at least two of the setsof seismic data, wherein the at least two sets of seismic data wererecorded at different times; applying the deconvolution operator in adeconvolution process to both of the at least two sets of data; andconducting further time-lapse processing to form a difference record. Insome cases, there are at least three sets of seismic data, wherein eachof the at least three set represents recordings from surveys conductedat different times, and wherein the designing of a deconvolutionoperator is dependent upon data from each of the sets of surveys. Insome of these cases, a further step is provided, comprising designing asecond deconvolution operator dependent upon data from the first and thethird sets of data, the second deconvolution operator being applied tothe first set of data and the third set of data. The first set of datamay have been taken before or after the second set of data with equaleffectiveness.

According to a more specific example embodiment, the design of the firstdeconvolution operator comprises: averaging of a power spectrum for afirst ensemble of traces from the first set; averaging of a powerspectrum for a second ensemble of traces from the second set; averagingof the power spectrum averages of the first and the second ensemble; anddesigning a deconvolution operator for the first and the second surveyfrom the averaging of the power spectrums of the first and the secondensemble.

In some circumstances, the above averaging for each subsequent surveywould be prohibitively expensive. Therefore, according to an evenfurther embodiment of the invention, multiple surveys are handled bystoring the averaging of the first survey, so that when a later surveyis made, the first survey's autocorrelations will not have to berecalculated. According to such an embodiment, the designing of thefirst deconvolution operator comprises: averaging of a power spectrumfor a first ensemble of traces from the first set; inverse transformingthe average into a time-domain representation of the average for thepower spectrum of the first ensemble; storing the time-domainrepresentation; averaging of a power spectrum for a second ensemble oftraces from the second set; inverse transforming the average into atime-domain representation of the average of the power spectrum for thesecond ensemble; averaging of the time-domain representations of theaverages for the power spectrum of the first and the second ensemble;and designing a deconvolution operator for the first and the secondsurvey from the average of the time-domain representations of theaverages for the power spectrum of the first and the second ensemble.

Referring now to FIG. 1A, a specific example will be discussed. FIG. 1Ashows a plot of a first raw stack (no DMO, migration, or other noisesuppression used, only spherical divergence and geometry correction) ofa survey taken at a first time is shown, without deconvolution. FIG. 1Bis a plot of the first record of FIG. 1A with deconvolution, wherein thedeconvolution operator was designed before the stack dependent upon thedata of the first raw stack. The specific parameters were:

Sample rate: 2 milliseconds

Opperator Length: 140 milliseconds

Prediction gap: Spiking, 2 milliseconds

Analysis Window: Offset dependent: for the near offsets, 300-5000 msec.,for the far offsets, 3600-5000 msec.

Added white noise: 0.5%

Decon apply window: 0-600 msec.

FIG. 2A is a plot of a second raw stack of a second survey taken of thesame geographic area as that of FIG. 1A, but at a different time,without deconvolution, and FIG. 2B is a plot of the second record ofFIG. 2A, with deconvolution, wherein the deconvolution operator wasdesigned before the stack dependent upon the data of the second rawstack. Here, the specific parameters were as above.

FIG. 3A is a plot of the difference in FIGS. 1A and 2A. FIG. 3B is aplot of the difference between FIGS. 1B and 2B. Ideally, there would benothing but noise or changes in structure showing in difference plot ofFIG. 3A, if the two surveys were otherwise identical and there had beenno change in the state of the reservoir. However, it can be seen thatunchanged structure appears in the difference plot 3B.

FIG. 1C is a plot of the first stack of FIG. 1A, with pre-stackdeconvolution, wherein the deconvolution operator was designed dependentupon the data of the first record and the second stack. FIG. 2C is aplot of the second stack of FIG. 2A, with pre-stack deconvolution,wherein the deconvolution operator was designed dependent upon the dataof the first record and the second record. The specific parameters were,for both FIGS. 1C and 2C, the same as for FIG. 1A.

FIG. 3C is a plot of the difference between the plots of FIGS. 1C and2C. Notice the lack of unchanged structure.

In designing the deconvolution filter used some embodiments, the first(i.e. “base”) dataset is processed, wherein the average of the traceautocorrelation functions within each ensemble (i.e. shot, receiver,CMP, etc.) is computed in the f-x domain, to estimate the powerspectrum:

i=0${{P\left( {i,k} \right)} = {\left( {1/N} \right){\sum\limits_{k = 1}^{N}{P\left( {i,k} \right)}}}},{i = 1},{\ldots \quad Q}$

where:

N is the number of traces in the x direction,

Q is the nyquist frequency index,

P(i, k) is the power spectrum estimation at frequency index i and tracenumber k, and

i is the time-lapse index.

Next, an inverse FFT is performed and the autocorrelation functiontraces for each ensemble are saved. The same power spectrum estimationis performed on a dataset taken at a different time (i.e. a “monitor”dataset) to compute:

i=0

P(i, k)

where i=1 is the first time lapse monitor survey.

Next, the all the time lapse autocorrelation average functions areaveraged to generate updated autocorrelation functions, according to thefollowing formula:${\overset{\_}{A}}_{BM} = {\sum\limits_{l = 1}^{M}\frac{P^{1}}{M}}$

i=1, . . . M where M is the number of time lapses.

The resulting autocorrelation traces are used to design convolutionfilters which are applied to all the time lapse surveys.

Other deconvolution filter design processes (e.g., single windowdeconvolution, time varying deconvolution, surface consistentdeconvolution, and common domain deconvolution (shot, receiver,midpoint, offset)) will occur to those of skill in the art, which aredependent upon both sets of data, which do not depart from the spirit ofthe present invention.

Further, like the deconvolution operation described above, there areother processing steps which benefit under one aspect of the presentinvention by designing the process parameters dependent upon the datafrom each set, rather than the date from one set. Examples include:crossequalization (discussed in more detail below), f-k filtering, τ-pand Radon domain filtering, static calculations, and multiple removal inFK, τ-p, and Radon domains.

There are always variations in actual acquisition, no matter how welloperations are planned or executed. And, in an ideal world theprocessing would also be identical, because the acquisition would bereplicated from base survey to monitor survey. But that is notrealistic. Thus, intersurvey acquisition variations occur with manysystems that are useful with the present invention (for example, fixedinstallation systems, OBS systems, and streamer systems). Thesevariations could be caused by feathering, cable drift, shooting azimuth,source variation, sea state conditions, etc., which can result inprocessing variations between surveys.

Crossequalization

According to one example embodiment of the invention, there are four“correction” elements the crossequalization function used. They aretiming corrections, RMS energy balancing, bandwidth normalizing, andphase matching. Each element essentially builds a transfer or impulseresponse function between the two surveys, or alternatively, between twotraces in each survey. The resultant, crossequalized trace (t_(XEQ)) iscalculated in the following manner:

t _(XEQ) =t*f(s _(corr′) rms _(corr′) m _(corr′) p _(corr′))  (1)

where t is the input trace, * connotes convolution, and s_(corr′),rms_(corr′), m_(corr′) and p_(corr′) are crossequalization elements (inimpulse response f) corresponding to time, amplitude, magnitude, andphase, respectively.

According to alternative embodiments of the invention, impulse responsefunctions are computed trace-by-trace, in a surface consistent manner,or globally from the base-line survey (base) and the monitor survey(repeat). According to the present example embodiment, a monitor surveyis transformed (crossequalized) to look like a base survey, and thecrossequalization impulse response functions are designed byhorizon-specific time gates that exclude the reservoir zone where changeis anticipated.

The effects of each component through simple models and example datafrom a time-lapse three-dimensional project will now be explained.

Crossequalization Elements

Time Corrections (s_(corr′))

Besides rough seas in a storm and current conditions that may causestreamer feathering, tidal and temperature changes are also in play.These changes may not be significant when recording, processing orinterpreting seismic data for a single 3D survey. However, when several3D surveys are involved for time-lapse monitoring, there are potentialintersurvey sea state variations.

Seasonal and shorter-period (storm related) variations in temperatureand salinity can also cause intersurvey differences between timereflectors. The degree of variation is a function of the severity of thechange between the temperature and salinity profiles and the depth ofthese variations in the water column. For example, surveys taken in theEast China sea show that over a 10 year period large scale, seasonalvariations are conspicuous, and they are sufficient to cause a temporalintersurvey variation of two milliseconds (TWT) between referencereflectors during seasonal extremes (if acquired in January and June,per se), where water depths are equal to or greater than 100 m.

As for tidal changes, an extreme example would be to shoot seismic datain the Bay of Fundy, between the Canadian Atlantic provinces of NewBrunswick and Nova Scotia. Here, the tidal changes are approximately 15m that is equivalent to 20 ms two-way time (assuming a typical watervelocity of 1500 m/s). If one survey was acquired at low tide andanother at high tide (assuming very short survey periods), seismicdifferencing of the datasets without a time correction would not yieldsatisfactory results.

In general, the errors increase with increasing delay, and broaderbandwidth wavelets obtain larger error differences at shorter delays.

RMS Energy Balancing (rms_(corr′))

In a perfect time-lapse seismic monitoring experiment, the same seismiccrew would acquire each monitor survey with the exact same equipment,under the exact (sea state) conditions that the original baseline surveywas acquired with. However, acquisition differences do occur, even underthe best controlled situations. For example, scaling the data toequivalent RMS levels is important, especially if it (the RMS energy)alters significantly between surveys.

In one example, from approximately 100 CMPs in a time-lapse 3D datasetacquired by streamer vessels approximately 18 months apart, acquisitiongeometry and instrumentation were nearly identical, with the exceptionthat one survey was acquired in the summer months and the other in thewinter. Analysis of the average amplitude of the seismic data withindiscrete 200 ms windows shows that there is amplitude bias between thebase survey and the monitor survey throughout the profile. Afterapplication of a smoothed, amplitude-scaling crossequalization elementwithin each 200 ms time gate, it was seen that the disparities betweenthe base and monitor surveys are smaller after crossequalization.Without question, variations such as these need to be corrected beforetime-lapse seismic difference sections/surveys can be interpreted.

Bandwidth Normalization (m_(corr′))

Equating bandwidths is also part of seismic differencing. If two seismicvolumes are to be subtracted, it would be best to do so with each volumehaving equivalent spectra. Crossequalization (if robust enough), shouldcorrect one spectrum to match another, but if the crossequalization isnot optimized, discrepancies appear. Larger errors occur for largemismatches in bandwidth. Tests indicate that residual reflector energy(i.e. energy left in the difference section/survey that, ideally, shouldbe removed) ranges from 7% to 40% of input amplitude for centerfrequency differences of 2.5 to 12.5 Hz, respectively (for Rickerwavelets with center frequencies incrementally separated by 2.5 Hzsubtracted from a baseline 30 Hz Ricker wavelet). A maximum centerfrequency spread of 12.5 Hz is reasonable for most vintage 3D surveys,but larger bandwidth difference occurs in some cases, and, in thosecases, the residual reflector energy is larger. Therefore, according toembodiments processing these surveys, crossequalization normalizes thebandwidth differences to reduce these potential seismic differencingerrors.

Phase Matching (p_(corr′))

Finally, errors are also associated with phase mismatches. For example,for a 30 Hz Ricker wavelet, phase-rotated by a constant 5, 10, 15, 30,60 and 90 degrees from zero-phase, residual reflector energy caused byphase differences can be as large as 20% of the input for phasemismatches as small as 15 degrees. Further, the human eye can havedifficulty seeing phase differences of 15 degrees or less, so it issignificant that the crossequalization operator correct for small phasemismatches as well as large ones. Phase difference errors are unaffectedby bandwidth differences. It has been found that residual reflectorerrors of 20% can occur for wavelets out of phase by 15 degrees.

Seismic Differencing

The ultimate goal in time-lapse seismic monitoring is to see changesthat infer interpore fluid movement or the absence of interpore fluidmovement between calendar dates. Currently this is obtained bysubtracting a monitor survey from the base survey after rectifying thedata through crossequalization. Differenced data with reduced amounts ofresidual reflector energy will have a higher probability of identifyingfluid movement than those with higher residual reflector energy. Tocharacterize the effectiveness of crossequalization, we show severaldifference plots with various component crossequalization applicationsusing a North Sea time-lapse experiment. Migrated inline stacks of thebase and monitor survey traversing a non-reservoir area are presented inFIGS. 4 and 5, respectively, and a prominent series of strong staticreflectors is visible between 2000 and 2600 ms. For this example adesign window of 1900 ms to 2700 ms is used. An 800 ms operator wascreated and applied from the design window, and a bandpass filter of3/8-35\55 Hz was applied to both datasets before crossequalization.Since these are static reflectors (i.e., they are reflection eventsrepresenting an unchanged portion of the geologic structure in or nearthe reservoir), subtraction of the base and monitor survey shouldideally result in very little residual reflector energy, since bothsurveys were acquired identically. However, the difference sectionobtained without any crossequalization, showed that subtraction did not.The residual reflector energy is significant. With these large amountsof residual reflector energy left in this section, identifying any fluidmovement (elsewhere) in reservoir portions of the volume would bedifficult.

Comparison of a difference section obtained without anycrossequalization with one in which only the amplitude spectrum wascrossequalized, one in which only the phase spectrum was crossequalized,and one in which both were crossequalized, shows that phase correctionsare a major component of the total crossequalization operation, whichreiterates the synthetic modeling results.

The above processing is performed on a massively parallel processorplatform (e.g. IBM SP2, Intel Paragon) and software compatible with suchhardware (e.g. PGS Tensor Inc.'s CUBE MANAGER™ software).

Suggested Readings

Two published case histories on time-lapse seismic monitoring areGreaves and Fulp (Geophysics, 1987), incorporated herein by reference,who monitor a fire-flood in west Texas, and Johnstad, Uden and Dunlop(First Break, 1993), also incorporated herein by reference, monitoringgas injection in the Oseberg Field located in the Norwegian sector ofthe North Sea. Another good paper involving a steam flood in the DuriField in Indonesia was presented by Lumley, et. al. at the 1995 SEGmeeting in Houston, which can be through the SEG.

The above description is given by way of example only and otherembodiments will occur to those of skill in the art without departingfrom the scope of the invention which is defined by the claims below andequivalents thereof.

What is claimed is:
 1. A method of comparing multiple seismic surveydata sets of a reservoir, wherein a first seismic survey data set istaken at a first time and a second seismic survey data set is taken at asecond time, to detect changes in the reservoir between the first timeand the second time, the method comprising: providing a first reflectionevent in the first seismic survey data set having a corresponding secondreflection event in the second seismic survey data set, wherein thefirst reflection event and the second reflection event represent anunchanged portion of geologic structure in or near the reservoir andwherein the first reflection event is represented by a first set ofevent parameters and the second reflection event is represented by asecond set of event parameters; providing an acceptance thresholddifference function between the first set of event parameters and thesecond set of event parameters; determining a crossequalization functionto apply to the second set of event parameters, the crossequalizationfunction being characterized in that, upon application of thecrossequalization function to the second set of event parameters,whereby a crossequalized set of event parameters is defined, and thedifference between the first set of event parameters and thecrossequalized set of event parameters is below the threshold differencefunction; applying the crossequalization function to a third reflectionevent, the third reflection event being related to the second data set,wherein a crossequalized third reflection event is defined, wherein thethird reflection event has a corresponding fourth reflection event inthe first data set, and wherein the third and fourth reflection eventsrepresent a changing portion of the reservoir; comparing thecrossequalized third reflection event to the fourth reflection event bysubtracting the crossequalized third reflection event from the fourthreflection event.
 2. A method as in claim 1 wherein said providing saidacceptance threshold difference function comprises: iterative selectionof event parameter modifications to the second set of event parameters,application of the event parameter modifications to the second set ofevent parameters, wherein a modified set of event parameters is defined,comparison of the modified set of event parameters to the first set ofevent parameters, wherein said iterative selection continues until aconvergence is reached, and wherein the acceptance threshold differencefunction comprises the modified set of event parameters at convergence.3. A method as in claim 2 wherein one of said event parameters compriseamplitude.
 4. A method as in claim 2 wherein one of said eventparameters comprise phase.
 5. A method as in claim 2 wherein one of saidevent parameters comprise bandwidth.
 6. A method as in claim 2 whereinone of said event parameters comprise time.
 7. A method as in claim 2wherein said event parameters comprise amplitude, phase, bandwidth, andtime.
 8. A method as in claim 1 wherein said determining acrossequalization function comprises: iterative selection of eventparameter modifications to the second set of event parameters,application of the event parameter modifications to the second set ofevent parameters, wherein a modified set of event parameters is defined,comparison of the modified set of event parameters to the first set ofevent parameters, and providing an acceptance threshold difference,wherein said iterative selection continues until a comparison resultfrom said comparison designates a difference between the first set ofevent parameters and the modified set of event parameters below theacceptance threshold difference.
 9. A method as in claim 1 wherein saidproviding an acceptance threshold difference function comprises:providing a windowed trace difference between a time window of a firsttrace from the first seismic survey data set and a time window of asecond trace from the second seismic survey data set, wherein the secondtrace includes reflection events corresponding to reflection events inthe first trace and wherein the time window of the second trace issubstantially the same as the time window of the first trace, andproviding a ratio of the windowed trace difference over the time windowof the first trace, choosing the acceptance threshold difference to beless than the ratio.
 10. A method as in claim 9 wherein the differencein bandwidth of the time window of the first trace and the time windowof the second trace is less than about 25%.
 11. A method as in claim 10wherein the difference in width of the time window of the second traceand the bandwidth of the time window of the first trace is less thanabout 10%.
 12. A method as in claim 1 wherein said providing anacceptance threshold difference function comprises: providing a windowedtrace difference between a time window of the square of a first tracefrom the first seismic survey data set and a time window of the squareof a second trace from the second seismic survey data set, wherein thesecond trace includes reflection events corresponding to reflectionevents in the first trace and wherein the time window of the secondtrace is substantially the same as the time window of the first trace,and providing a ratio of the windowed trace difference over the timewindow of the square of the first trace, choosing the acceptancethreshold difference to be less than the ratio.
 13. A method as in claim12 wherein the time window has a length of at least about two reflectionevents.
 14. A method as in claim 13 wherein the time window has a lengthof at least about five reflection events.
 15. A method as in claim 1wherein said applying the crossequalization function to a thirdreflection event in the second data set comprises convolution betweenthe crossequalization function and the third reflection event in thesecond data set.
 16. A method as in claim 1 wherein said first data setcomprises a trace from a seismic receiver.
 17. A method as in claim 1wherein said first data set and said second data set comprise a summedset of traces from a set of seismic receivers.
 18. A method as in claim1 wherein said first data set and said second data set comprise a summedset of traces from a set of borehole receivers.
 19. A method as in claim1 wherein said first data set and said second data set comprise prestackdata.
 20. A method as in claim 19 wherein said prestack data comprisesCMP data.
 21. A method as in claim 19 wherein said prestack datacomprises shot data.
 22. A method as in claim 19 wherein said prestackdata comprises migrated data.
 23. A method as in claim 1 wherein saidfirst data set and said second data set are subjected to equivalentprestack processes.
 24. A method as in claim 1 wherein the first dataset is subjected to the same designature process as the second data set.25. A method as in claim 1 wherein the first data set is subjected tothe same deconvolution process as the second data set.
 26. A method asin claim 1 wherein the first data set is subjected to the same noiseattenuation processing steps as second data set.
 27. A method as inclaim 1 wherein the first data set is subjected to the same multipleattenuation processing as the second data set.
 28. A method as in claim1 wherein the same DMO operator is used on first data set as on thesecond data set.
 29. A method as in claim 1 wherein migration on thefirst data set is conducted with the same velocity field as migration onthe second data set.
 30. A method as in claim 1 wherein migration on thefirst data set is conducted with the same migration operator asmigration on the second data set.
 31. A method as in claim 1 whereinfiltering on the first data set is conducted with the same filter asfiltering on the second data set.
 32. A method of deconvolution ofmultiple sets of seismic data from the same geographic area is provided,the method comprising: designing of a deconvolution operator dependentupon data from at least two of the sets of seismic data, wherein the atleast two sets of seismic data were recorded at different times;applying the deconvolution operator in a deconvolution process to bothof the at least two sets of data; and conducting further time-lapseprocessing to form a difference record.
 33. A method as in claim 32,wherein the at least two of the sets of data comprises at least threesets of seismic data, wherein each of the at least three sets representsrecordings from surveys conducted at different times, and wherein thedesigning of a deconvolution operator is dependent upon data from eachof the sets of surveys.
 34. A method as in claim 32, wherein the atleast two of the sets of data comprises at least three sets of seismicdata, wherein each of the at least three sets represents recordings fromsurveys conducted at different times, and wherein the designingcomprises designing a first deconvolution operator dependent upon datafrom a first set of data and a second set of data, the firstdeconvolution operator being applied to the first set of data and thesecond set of data, and further comprising: designing a seconddeconvolution operator dependent upon data from the first and a thirdset of data, the second deconvolution operator being applied to thefirst set of data and the third set of data.
 35. A method as in claim34, wherein the first set of data represents recordings from a surveytaken before a survey represented by the second or third sets of data.36. A method as in claim 34, wherein the designing of the firstdeconvolution operator comprises: averaging of a power spectrum for afirst ensemble of traces from the first set; inverse transforming theaverage into a time-domain representation of the average for the powerspectrum of the first ensemble; storing the time-domain representation;averaging of a power spectrum for a second ensemble of traces from thesecond set; inverse transforming the average into a time-domainrepresentation of the average of the power spectrum for the secondensemble; averaging of the time-domain representations of the averagesfor the power spectrum of the first and the second ensemble; anddesigning a deconvolution operator for the first and the second surveyfrom the average of the time-domain representations of the averages forthe power spectrum of the first and the second ensemble.
 37. A method asin claim 32, wherein the designing of the a deconvolution operatorcomprises: averaging of a power spectrum for a first ensemble of tracesfrom the first set; averaging of a power spectrum for a second ensembleof traces from the second set; averaging of the power spectrum averagesof the first and the second ensemble; and designing a deconvolutionoperator for the first and the second survey from the averaging of thepower spectrums of the first and the second ensemble.
 38. A method forperforming time-lapse seismic survey, signal processing, having at leasta first survey and a second survey, wherein each survey has a set ofseismic data, the method comprising: performing a set of processingsteps on the first survey; performing the set of processing steps on thesecond survey; wherein the set of processing steps is dependent upon aset of seismic signal parameters; choosing at least one of the set ofparameters by a selection process dependent upon data from both thefirst survey and the second survey; and applying the at least one of theset of parameters in the at least one of the set of processing steps toboth the first survey and the second survey.
 39. A method as in claim 38wherein the at least one of the sets of parameters comprises adeconvolution operator.
 40. A method as in claim 39 further comprising:designing of a deconvolution operator dependent upon data from at leasttwo of the sets of seismic data, wherein the at least two sets ofseismic data were recorded at different times; applying thedeconvolution operator in a deconvolution process to both of the atleast two sets of seismic data; and forming a difference recorddependent upon the applying the deconvolution operator in thedeconvolution process to both of the at least two sets of seismic data.41. A method as in claim 40, wherein the at least two of the sets ofseismic data comprises at least three sets of seismic data, wherein eachof the at least three sets represents recordings from surveys conductedat different times, and wherein the designing of a deconvolutionoperator is dependent upon data from each of the at least three sets ofseismic data.
 42. A method as in claim 40, wherein the at least two ofthe sets of seismic data comprises at least three sets of seismic data,wherein each of the at least three sets represents recordings fromsurveys conducted at different times, and wherein the designingcomprises designing a first deconvolution operator dependent upon datafrom a first set of data and a second set of data, the firstdeconvolution operator being applied to the first set of data and thesecond set of data, and further comprising: designing a seconddeconvolution operator dependent upon data from the first set of dataand a third set of data, the second deconvolution operator being appliedto the first set of data and the third set of data.
 43. A method as inclaim 42, wherein the first set of data represents recordings from asurvey taken before a survey represented by the second or third sets ofdata.
 44. A method as in claim 42, wherein the designing of the firstdeconvolution operator comprises: averaging of a power spectrum for afirst ensemble of traces from the first set; inverse transforming theaverage into a time-domain representation of the average for the powerspectrum of the first ensemble; storing the time-domain representation;averaging of a power spectrum for a second ensemble of traces from thesecond set; inverse transforming the average into a time-domainrepresentation of the average of the power spectrum for the secondensemble; averaging of the time-domain representations of the averagesfor the power spectrum of the first and the second ensemble; anddesigning a deconvolution operator for the first and the second surveyfrom the average of the time-domain representations of the averages forthe power spectrum of the first and the second ensemble.
 45. A method asin claim 40, wherein the designing of the deconvolution operatorcomprises: averaging of a power spectrum for a first ensemble of tracesfrom the set of data for the first survey; averaging of a power spectrumfor a second ensemble of traces from the set of data for the secondsurvey; averaging of the power spectrum averages of the first and thesecond ensemble; and designing a deconvolution operator for the firstand the second survey from the averaging of the power spectrums of thefirst and the second ensemble.